Модели стационарных временных рядов

В этой лекции через \(y_t\) будет обозначать наблюдаемый временной ряд, а \(\epsilon_t\) ненаблюдаемый процесс белого шума. Общим линеным процессом типа скользящего среднего называется процесс вида \(y_t=\epsilon_t+\Psi_1\epsilon_{t-1}+\Psi_2\epsilon_{t-2}+....+\Psi_n\epsilon_{t-n}+....\) Введем линейный оператор сдвига назад \(B: By_t=y_{t-1}\)

По индукции \(B^k: B^ky_t=y_{t-k}\)

Используя оператор сдвига, определение линеным процесса типа скользящего среднего будем записывать в операторном виде \(y_t=(1+\Psi_1B+\Psi_2B^2+.....)\epsilon_t= \Psi(B)\epsilon_t\)

Будем предполагать, что

\(\sum_{i=1}^{\infty}\Psi_i < \infty\)

Важный нетривиальный пример такого процесса это процесс, где все \(\Psi_i=\phi^i\), при \(-1<\phi<1\)

В этом случае

\(y_t= (1+\phi B+\phi^2B^2+...)\epsilon_t=\epsilon_t+\phi\epsilon_{t-1}+\phi^2\epsilon_{t-2}+...\)

Так как \(\epsilon_t\) - процесс белого шума \(E[\epsilon_t]=0\) поэтому

\(E[y_t]=E[\epsilon_t+\phi\epsilon_{t-1}+\phi^2\epsilon_{t-2}+...]= 0\)

Пусть \(D[\epsilon_t]=\sigma_{\epsilon}^2\) тогда

\(D[y_t]=D[\epsilon_t+\phi\epsilon_{t-1}+\phi^2\epsilon_{t-2}+...]= \sigma_{\epsilon}^2+\phi^2\sigma_{\epsilon}^2+\phi^4\sigma_{\epsilon}^2+... = \sigma_{\epsilon}^2(1+\phi^2+\phi^4+...)= \frac{\sigma_{\epsilon}^2}{1-\phi^2}\)

Вычислим автоковариационную и автокорреляционную функции процесса

\(c_1=cov(y_t,y_{t-1})= cov(\epsilon_t+\phi\epsilon_{t-1}+\phi^2\epsilon_{t-2}+...,\epsilon_{t-1}+\phi\epsilon_{t-2}+\phi^2\epsilon_{t-3}+...)= \phi\sigma_{\epsilon}^2(1+\phi^2+\phi^4+...)=\frac{\phi\sigma_{\epsilon}^2}{1-\phi^2}\)

\(corr(y_t,y_{t-1})= \frac{\phi\sigma_{\epsilon}^2}{1-\phi^2}/(\frac{\sigma_{\epsilon}^2}{1-\phi^2})= \phi\)

Аналогично вычислим

\(c_k=cov(y_t,y_{t-k})=\frac{\phi^k\sigma_{\epsilon}^2}{1-\phi^2}\)

\(\gamma_k=corr(y_t,y_{t-k})=\phi^k\)

Для произвольного линейного процесса скользящего среднего \(y_t=(1+\Psi_1B+\Psi_2B^2+.....)\epsilon_t= \Psi(B)\epsilon_t\) Нетрудно вычислить, что

\(E[y_t=0]\),

\(c_k=cov(y_t,y_{t-k})=\sigma_{\epsilon}^2\sum_{i=0}^{\infty}\Psi_i\Psi_{i+k}\)

Важно отметить, что ковариационная функция не зависит от \(t\) , зависит только от задержки \(k\) cледовательно это стационарный процесс с нулевым среднем. Для получения процесса с произвольным среднем \(\mu\) просто добавим \(\mu\) к правой части соотношения \(y_t=(1+\Psi_1B+\Psi_2B^2+.....)\epsilon_t+\mu= \Psi(B)\epsilon_t+\mu\)

В случае конечного числа \(\Psi_k\) отличных от нуля модель записывают в виде

\(y_t=\mu+\epsilon_t-\theta_1\epsilon_1- ...-\theta_q\epsilon_q\)

называют моделью скользящего среднего порядка \(q\) и обозначают \(MA(q)\)

Модель \(MA(1)\)

Подробно рассмотрим модель скользящего среднего первого порядка с нулевым среднем \(y_t=\epsilon_t-\theta\epsilon_{t-1}\)

Cогласно проведенным выше вычислениям сразу запишем, что автоковариационная функция процесса \(y_t\)

\(c_0=D[y_t]= \sigma_{\epsilon}^2(1+\theta^2)\)

\(c_1=cov(y_t,y_{t-1})=-\sigma_{\epsilon}^2\theta\)

\(c_k = 0\) при \(|k|>1\)

автокорреляционная функция процесса \(y_t\)

\(\rho_0=1\)

\(\rho_1= -\theta/(1+\theta^2)\)

\(\rho_k= 0\) при \(|k|>1\)

Поведение процесса, выборочной и модельной (True ACF) автокорреляционной функции при различных параметрах \(\theta\) можно изучить в прилагаемом ниже фрейме. Понятие частной автокорреляционной функции (PACF) и спектра в курсе дано будет позднее

Пример моделирования MA(1) на R при \(\theta=0.9\)

library(fArma)

## Loading required package: timeDate

## Loading required package: timeSeries

## Loading required package: fBasics

##

## Rmetrics Package fBasics

## Analysing Markets and calculating Basic Statistics

## Copyright (C) 2005-2014 Rmetrics Association Zurich

## Educational Software for Financial Engineering and Computational Science

## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.

## https://www.rmetrics.org --- Mail to: info@rmetrics.org

library(TSA)

## Loading required package: leaps

## Loading required package: locfit

## locfit 1.5-9.1    2013-03-22

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-9. For overview type 'help("mgcv-package")'.

## Loading required package: tseries

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:timeDate':
## 
##     kurtosis, skewness

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

theta <- 0.9
modelMA <- list(ma= theta)
ma1<-armaSim(modelMA,n=200)
matplot(ma1,main= "Simulated MA(1) process",type = "l",lwd=3,col="blue")

Для того чтобы визуально увидеть зависимость ряда от прошлого полезно строить график \(y_t\) от \(y_{t-1}\) и увидеть положительную корреляцию между ними

plot(ma1[,1,1],zlag(ma1[,1,1]), ylab = expression(y[t]),xlab = expression(y[t-1]),type='p',col = 'blue')

Корреляции между \(y_t\) от \(y_{t-2}\) уже не будет.

plot(ma1[,1,1],zlag(ma1[,1,1],2), ylab = expression(y[t]),xlab = expression(y[t-2]),type='p',col = 'blue')

Процесс MA(q)

Теперь рассмотрим процесс скользящего среднего произвольного порядка \(MA(q)\)

\(y_t=\epsilon_t-\theta_1\epsilon_{t-1}-...-\theta_q\epsilon_{t-q}\)

Для процесса \(y_t\) дисперсия

\(c_0=D[y_t]=\sigma_{\epsilon}^2(1+\theta_1^2+....+\theta_q^2)\)

автокорреляционная функция

\(\rho_k= \frac{-\theta_{k}+\theta_{k-1}\theta_{1}-\theta_{k-2}\theta_{2}+\theta_{k-q}\theta_{q}}{(1+\theta_1^2+....+\theta_q^2)}\) при \(k=1,...,q\)

\(\rho_k=0\) при \(k>q\)

Поведение процесса, выборочной и модельной (True ACF) автокорреляционной функции при различных параметрах \(\theta_1,\theta_2,\theta_3,\) можно изучить в прилагаемом ниже фрейме на примере процесса \(MA(3)\).

Процессы авторегрессии

Название происходит от регрессии самого на себя. Процесс вида

\(y_t= \phi_1y_{t-1}+\phi_2y_{t-2}+...+\phi_py_{t-p}+\epsilon_t\)

назвается процессом авторегрессии порядка \(p\) обозначается \(AR(p)\). Запись с использованием оператора сдига назад \(B\) и в операторной форме

\((1 - \phi_1B+\phi_2B^2-...- \phi_pB^p) y_t=\phi(B)y_t=\epsilon_t\)

Процесс авторегрессии 1 порядка AR(1)

Пусть стационарный процесс c нулевым среднем представим в виде

\(y_t= =\phi y_{t-1}+\epsilon_t\)

или, используя оператор сдвига

\(\phi(B)x_t = \epsilon_t\)

где \(\phi(B)=1-\phi B\)

Вычислим дисперсию \(c_0=D[y_t]\)

\(c_0=\phi^2c_0+\sigma_{\epsilon}^2\)

Умножая выражение для процесса на \(y_{t-k}\) и ,беря математическое ожидание от обеих частей получим выражение для ACF процесса

\(E[y_ty_{t-k}]= \phi E[y_{t-1}y_{t-k}]+E[\epsilon_ty_{t-k}]\)

\(c_k= \phi c_{k-1}\)

Откуда

\(c_k=\phi^k \frac{\sigma_{\epsilon}^2}{1-\phi^2}\)

Автокорреляционная функция

\(\rho_k=\frac{c_k}{c_0}=\phi^k\)

Предположим, чьл \(|\phi |< 1\) в этом случае \(\rho_k->0\) при \(k->\infty\)

Поведение процесса, выборочной и модельной (True ACF) автокорреляционной функции при различных параметрах \(\phi\) можно изучить в прилагаемом ниже фрейме на примере процесса \(AR(1)\).

Пример моделирования AR(1) на R при \(\phi=0.9\)

phi <- 0.9
model <- list(ar= phi)
ar1<-armaSim(model,n=200)
matplot(ar1,main= "Simulated AR(1) process",type = "l",lwd=3,col="blue")

Процесс AR(1) стационарен тогда и только тогда, \(|\phi|<1\), что эквиваленетно тому, что корень уравнения \(\phi(B) = 0\) по модулю больше 1.

plot(ar1[,1,1],zlag(ar1[,1,1]), ylab = expression(y[t]),xlab = expression(y[t-1]),type='p',col = 'blue')

Остается корреляция между \(y_t\) от \(y_{t-2}\)

plot(ar1[,1,1],zlag(ar1[,1,1],2), ylab = expression(y[t]),xlab = expression(y[t-2]),type='p',col = 'blue')

Процесс авторегрессии 2 порядка. AR(2)

Пусть процесс c нулевым среднем представим в виде

\(y_t=\phi_1 y_{t-1}+\phi_2 y_{t-2}+ \epsilon_t\)

или, используя оператор сдвига

\(\phi (B)x_t = \epsilon_t\)

где \(\phi(B)=1-\phi_1 B- \phi_2 B^2\)

Можно покачать, что для стационарности AR(2) процесса необходимо и достаточно, чтобы корни характеристического уравнения \(\phi(B) = 0\) или \(1-\phi_1 B- \phi_2 B^2=0\) по модулю были больше 1. Корни характеристического уравнения легко найти

\(B_{1,2}=\frac{\phi_1+-\sqrt{\phi_1^2+4\phi_2}}{-2\phi_2}\)

Абсолютные значения корней должны быть больше 1. Отсюда можно получить условия на \(\phi_1,\phi_2\) а именно

\(\phi_1+\phi_2<1\),\(\phi_2-\phi_1<1\) и \(|\phi_2|<1\)

Вычислим автоковариационную и автокорреляционную функции процесса AR(2) Умножая предствление процесса на \(y_{t-k}\) и, беря математическое ожидание, получим

\(c_k= \phi_1c_{k-1}+\phi_2c_{k-2}\) \(:k=1,2,....\)

автокорреляционная функция

\(\rho_k= \phi_1\rho_{k-1}+\phi_2\rho_{k-2}\) \(:k=1,2,....\)

Последнее уравнение носит название уравнение Юла-Уокера. Его обычно решают численно, но существует и другой способ нахождения \(\rho_k\) обозначим через \(G_1,G_2\) корни характеристического уравнения \(1-\phi_1 B- \phi_2 B^2=0\)

\(G_1= \frac{\phi_1+\sqrt{\phi_1^2+4\phi_2}}{2}\)

\(G_2= \frac{\phi_1-\sqrt{\phi_1^2+4\phi_2}}{2}\)

Если \(G_1\ne G_2\) то можно показать,что

\(\rho_k=\frac{(1-G_2^2)G_1^{k+1}- (1-G_1^2)G_2^{k+1}}{(G_1-G_2)(1+G_1G_2)}\)

Eсли корни комплексны, тогда

\(\rho_k=R^k\frac{sin(\Theta k+\Phi)}{sin\Phi}\)

где \(R =\sqrt{-\phi_2}\), \(cos\Theta=\phi_1/(2\sqrt{-\phi_2})\) и \(tan\Phi=(1-\phi_2)/(1+\phi_2)\)

Поведение процесса, выборочной и модельной (True ACF) автокорреляционной функции при различных параметрах \(\phi_1,\phi_2\) можно изучить в прилагаемом ниже фрейме на примере процесса \(AR(2)\).

Пример моделирования AR(2) на R при \(\phi_1=0.4,\phi_2=-0.5\)

phi <- c(0.4,-0.5)
model <- list(ar= phi)
ar2<-armaSim(model,n=200)
matplot(ar2,main= "Simulated AR(2) process",type = "l",lwd=3,col="blue")

Дисперсия процесса AR(2)

Чтобы вычислить дисперсию процесса AR(2) возьмем оператор дисперсии от обеих частей \(y_t=\phi_1 y_{t-1}+\phi_2 y_{t-2}+ \epsilon_t\)

\(c_0=D[y_t]=(\phi_1^2+phi_2^2)+2\phi_1\phi2+\sigma_{\epsilon}^2\) \((*)\)

У нас есть уранение \(c_k= \phi_1c_{k-1}+\phi_2c_{k-2}\) \(:k=1,2,....\) при \(k=1\) оно примет вид \(c_1= \phi_1c_{0}+\phi_2c_{1}\) Уравненике \((*)\) и последнее уравнение нам даст

\(c_0=\frac{1-\phi_1}{1+\phi_2}\frac{\sigma_{\epsilon}^2}{(1-\phi_2)^2-\phi_1^2}\)

Процесс авторегрессии произвольного порядка. AR(p)

Процесс вида

\(y_t=\phi_1 y_{t-1}+\phi_2 y_{t-2}+...-+\phi_p y_{t-p}+ \epsilon_t\) \((1)\)

называется процессом авторегрессии порядка \(p\)

В операторной форме

\(\phi (B)y_t = \epsilon_t\)

\((1-\phi_1 B+...+\phi_p B^p)y_t = \epsilon_t\)

характеристическое уравнение

\(1-\phi_1 B+...+\phi_p B^p= 0\)

Процесс авторегрессии порядка \(p\) стационарен тогда и только тогда корни характеристического уравнения по модулю больше единицы, или корни лежат вне единичной окружности \(|z|= 1\)

Пример

\(\phi_1=0.6,\phi_2=-0.3,\phi_3= 0.2\)

phi1<-0.6
phi2<--0.3
phi3 <- 0.2
par(mfrow = c(1, 1), cex = 0.9)
armaRoots(c(phi1,phi2,phi3),lwd = 8, n.plot = 400, digits = 8)

##            re        im     dist
## 1  1.59042489  0.000000 1.590425
## 2 -0.04521244  1.772504 1.773080
## 3 -0.04521244 -1.772504 1.773080

Умножая (1) на \(y_{t-k}\), и ,разделив на \(c_0=D[y_t]\), берем математическое ожидание от обоих частей полученного соотношения, получим

\(\rho_k=\phi_1\rho_{k-1}+...+\phi_p \rho_{k-p}\)

Полагая последовательно \(k=1,..,p\), получим систему линейных алгебраических уравнений

\(\rho_1=\phi_1+...+\phi_p \rho_{p-1}\)

\(\rho_2=\phi_1\rho_1+...+\phi_p \rho_{p-2}\)

\(...\)

\(\rho_p=\phi_1\rho_{p-1}+...+\phi_p\)

Эта система уравнений носит название Юла-Уокера. Зная набор чисел \(\phi_1,\phi_2,...,\phi_p\) можно из системы вычислить \(\rho_1,\rho_2,...,\rho_p\) и наоборот.

С уравнением Юла Уокера связано появление функции частной автокорреляции (PACF). Последовательно будем решать их для \(k=1,2,...,p\) При \(k=1\) оно имет вид \[\rho_1=\phi_1\] Решение его обозначим \(\phi_{1,1}\).

Для \(k=2\) \[\rho_1=\phi_1+\phi_2\rho_{1}\]

\[\rho_2=\phi_1\rho_1+\phi_2\] Решение его обозначим \(\phi_{2,1},\phi_{2,2}\). И так далее повышаем порядок уравнения и решения для порядка \(k\) его обозначаем \(\phi_{k,1},\phi_{k,2},..., \phi_{k,k}\). Функция от \(k\), поученная из решений уравнения при старших степенях уравнений \[\phi_{1,1},\phi_{2,2},...,\phi_{p,p},...,\phi_{k,k}\] получила название частной автокорреляционной функции (PACF). Для AR(p) она по построению имеет отличные от нуля значения до \[\phi_{1,1},\phi_{2,2},...,\phi_{p,p}\] далее при \(k>p\) должны следовать нули так как соответствующие \(\phi_k\) в модели AR(p) имеются в наличии только до порядка \(p\). На практике конечно нули мы не увидим однако, дисперсия \[D[\phi_{k,k}]\approx\frac{1}{n}\] попадание значений PACF в полосу \(\pm\frac{1}{\sqrt{n}}\) является сигналом к идентификации порядка модели \(p\).

Выражение для дисперсии процесса получим из (1), умножив его на \(y_t\), и , беря математическое ожидание, получим

\(D[y_t]= c_0= \frac{\sigma_{\epsilon}^2}{1-\phi_1\rho_1-... -\phi_p\rho_p}\)

Поведение процесса, выборочной и модельной (True ACF,PACF) автокорреляционной функции и частной автокорреляционной при различных параметрах \(\phi_1,\phi_2,\phi_3\) можно изучить в прилагаемом ниже фрейме на примере процесса \(AR(3)\).

Идентификация моделей AR(p) и МА(q)

PACF и ACF служат прекрасным средством для идентификации модели и их порядков на практике. Сведем наши знания об этих функциях в следующую таблицу

  par(mfrow = c(2, 2), cex = 1.0,col = "blue",lwd = 3)
  phi1 = 0.6
  phi2 = -0.3
  phi3 = 0.3
  modelAR3 <- list(ar=c(phi1,phi2,phi3))
  armaTrueacf(modelAR3,lag.max = 15,type = "both", doplot = TRUE) 
  theta1 = 0.5
  theta2 = 0.3
  theta3 = 0.4
  modelMA3 <- list(ma=c(theta1,theta2,theta3))
  armaTrueacf(modelMA3,lag.max = 15,type = "both", doplot = TRUE)

Таблица почему-то не скомпилилась в html файл. Поэтому,использую операцию copy-paste перенесем R код в RStudio и посмотрим таблицу там.

Потренировать умение правильно идентификацировать заведомо известную модель AR(3) или MA(3) по выборочной ACF,PACF поможет следующий фрейм.

Смешаный процесс ARMA(p,q)

Процесс с нулевым среднем вида

\(y_t= \phi_1y_{t-1}+ ...+ \phi_py_{t-p}+\epsilon_t- \theta_1\epsilon_{t-1}-..-\theta_q \epsilon_{t-q}\)

называется cмешаным процессом авторегрессии-скользящего среднего порядка \(p,q\), обозначается \(ARMA(p,q)\) В операторной записи

\((1-\phi_1B-...-\phi_pB^P)y_t= (1-\theta_1B-...-\theta_pB^q)\epsilon_t\)

\(\phi_p(B)y_t= \theta_q(B)\epsilon_t\)

ARMA(1,1) модель

Начнем с простейшей модели ARMA(1,1).

\(y_t= \phi_1y_{t-1}+\epsilon_t- \theta_1\epsilon_{t-1} (2)\)

Так как

\(E\epsilon_t y_t= E[\epsilon_t\phi_1y_{t-1}]+E[\epsilon_t\epsilon_t]- E[\theta_1\epsilon_t\epsilon_{t-1}]= \sigma_{\epsilon}^2\)

\(E\epsilon_{t-1} y_t= E[\epsilon_{t-1}\phi_1y_{t-1}]+E[\epsilon_{t-1}\epsilon_t]- E[\theta_1\epsilon_{t-1}\epsilon_{t-1}]= \phi_1\sigma_{\epsilon}^2-\theta_1\sigma_{\epsilon}^2= \sigma_{\epsilon}^2(\phi_1-\theta_1)\)

Если умножить на \(y_{t-k}\) и взять математическое ожидание , то прдем к следуюшей системе соотношений \[c_0=\phi_1c_1+[1-\theta_1(\phi_1-\theta_1)]\sigma_{\epsilon}^2\] \[c_1= \phi_1c_0-\theta_1\sigma_{\epsilon}^2\] \[c_k= \phi_1c_{k-1}\]

Разрешая это систему получим, что \[\rho_k= \frac{(1-\theta_1\phi_1)(\phi_1-\theta_1)}{1-2\theta_1\phi_1-\theta_1^2}/phi_1^k\] для \(k>0\).

Таким образом ACF убывает с ростом \(k\) также как и модель AR(1)

ARMA(p,q) модель

Запишем модель в операторном виде. \[\phi_p(B)y_t=\theta_q(B)\epsilon_t\] Без доказательства (громоздкие формулы) примем следующий факт, что в ACF первый \(q\) значений ведут себя произвольно,а начиная с задержки \(q+1\) идет экспоненциальное затухание аналогичное модели AR(p).

Стационарность модели зависит от характеристического многочлена \(\phi_p(B)=0\) в том смысле,что корни его должны лежать вне единичного круга \(|B| = 1|\).

Обратимость определяется многочленом \(\theta_q(B)=0\).Корни его соответственно должны лежать вне единичного круга \(|B| = 1|\).

Идентификация ARMA(p,q) моделей. EACF(Extended ACF)

Идентификация смешаных моделей только по ACF и PACF затруднительна и может привести к ошибке. Для примера приведем модель ARMA(1,1). А библилотеке TSA (версия 3.2.3) моделируется ARMA(1,1) процесс вид которого

  library(TSA)
  data("ma1.1.s")
  matplot(ma1.1.s, lwd = 2,col = "blue",type ="l")

Его ACF и PACF выглядят так

  acf(ma1.1.s)

  pacf(ma1.1.s)

Что позволяет выдвинуть гипотезу о модели AR(1).
Однако анализ остатков (в следующей лекции о нем) показывает, что в остатках не белый шум. Для облегчения идентификации смешаных моделей Tsau и Tiao в 1984 году предложили метод идентификации, получивший название Extended ACF. Суть его в одновременной проверке значимости оценок ACF,PACF. По результатам которого выводится таблица крестиков и ноликов. По левому верхнему углу устойчивой границы между которыми можно сделать вывод о предполагаемых порядках модели. Конечно только анализ остатков после оценки модели поможет сделать окончательный вывод

  eacf(ma1.1.s)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o x o o o o o o  o  o  x 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x x x x o o o o o o o  o  o  o 
## 4 x x o x x o o o o o o  o  o  o 
## 5 x x o o x x o o o o o  o  o  o 
## 6 o o x x x x o o o o o  o  o  o 
## 7 o o x o o x o o o o o  o  o  o

Для практического ознакомления со смешаными моделями предлагаю следующий фрейм.

Lecture 3

Doynikov A

December 11, 2015